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Abstract 

We review the diffuse scattering and the loss coefficient in ultracold neutron reflection from 
slightly rough surfaces, report a surprising reduction in loss coefficient due to roughness, and 
discuss the possibility of transition from quantum treatment to ray optics. The results are used in 
a computer simulation of neutron storage in a recent neutron lifetime experiment that reported a 
large discrepancy of neutron lifetime with the current particle data value. Our partial re-analysis 
suggests the possibility of systematic effects that were not included in this publication. 
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I. INTRODUCTION 



The scattering and absorption of cold and ultracold neutrons (UCN) at slightly rough 
surfaces resembles the roughness problem in light and X-ray optics [1]. Surface roughness 
induces losses in neutron guide tubes and affects the behavior of a UCN gas in a trap. 
We investigate especially the relation between a realistic account of roughness in computer 
simulations of UCN storage and a reliable data interpretation with respect to the neutron 
lifetime. The issue has gained added importance since a new UCN storage based lifetime 
value six standard deviation away from the world average value [2] was published in 2005 



Most theoretical analyses of diffuse cold and ultracold neutron scattering at surfaces 
with small roughness (e.g., [4-7]) are based on first-order perturbation theory. A possible 
extension to macroscopic ray optics was presented in [4]. Carrying the analysis to second 
order perturbation theory, Ignatovich [7] derived the roughness effect on the loss coefficient 
for reflection at a slightly absorbing wall material. The analysis was complex, and the 
result has, apparently, not been numerically evaluated and applied so far. We used a more 
direct approach, verified and quantified the result of [7], and in the process obtained a 
derivation of the " Debye- Waller factor" (describing the attenuation of the specular beam due 
to roughness), that is consistent with standard expressions but requires fewer assumptions 
about the roughness characteristics. This appears important since in practical cases very 
little is known about the roughness parameters of a given surface. Finally, we address the 
question to what extent a transition to a macroscopic picture is possible. 

II. CORRELATION FUNCTIONS DESCRIBING A ROUGH SURFACE 

Small irregular deviations of a slightly rough surface from the plane geometry are usually 
described by a height-height correlation 



where £(p) is the random elevation at point p = (x, y) of the plane surface above its average 
z — 0. A is the illuminated surface area and 6 is the displacement vector between two points. 
Among common models [4-6] we will concentrate on a Gaussian correlation for solid surfaces 



[3]. 
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with mean-square roughness b 2 = (£ 2 ): 



f G (6) = f G (6)=b 2 eMS 2 /(2w 2 )] (2) 

and a L K -modeV [5] for liquids and possibly glasses retaining characteristics of a liquid 
below the glass transition: 

tf„[V(W2)] 

The latter form contains the modified Bessel function K and has been proposed in [6] as 
a modification of a liquid model with logarithmic short-range divergence [5] (which would 
imply a divergence of g(S), see below) to account for smoothing due to surface tension. 
Smoothing is achieved by applying a short-range cutoff, So, in addition to the long-range 
cutoff w that is used for both models. In [4, 6] we have emphasized the importance of the 
slope-slope correlation for the surface gradient x — V£(p) : 

g(S) = g(S) = Jim - [ X (p) ■ X (p + 6) d 2 6 = -V 2 f(5) (4) 

A^oo A J A 

which is determined by f(S) through the Laplace operator V 2 . This can be verified using 
Gauss's divergence theorem. For the Gaussian model: 

g G (S) = a 2 G (l - exp[-5 2 /(2 W 2 )] (5) 

with mean-square slope Oq = (x 2 ) = g(0) = 2b 2 /iv 2 . For the 'K -model' [6] 



2 P{[(S 2 + S 2 )/(2u?p 2 } 



P[So/{wV2)\ 



with 

P(v) = s Mfl Mv) _ KoW) 



v 

and 

a 



2 _ I'' I^C) 



K w 2 tK (t) 

To facilitate comparison with the Gaussian model we choose Sq such that for given b and 
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w the mean-square slopes become identical: Oq = a 2 K = a 2 . This requires t = 5 /(w\/2) = 
0.7709 which is the solution of K 2 (t)/K (t) = 5. 

Fairly smooth surfaces with small a 2 are best characterized by g(S). For surfaces with 
tips and sharp edges the curvature-curvature correlation function 

h(S) = h{5) = Jim I jT k{ P )k{ P + 6) d 2 8 = iv 2 V 2 /W (7) 



may not be negligible, either. For a 2 « 1 the mean surface curvature is given by k — 
\ (§i? + f^) = an< ^ re P resen ts the mean curvature in any two orthogonal in-plane 

directions x and y. The operations in (7) can be performed using the properties of Bessel 
functions. The results are given in Appendix A. The mean values are 

4 = ft G ( ) = 2b 2 /w A = a 2 /v? 

and 

K 2 K = h K (0) = 2.103a 2 /^ 2 

for t = 0.7709, as before. 

We add a note related to calculations of 'Debye- Waller factors' describing the decrease of 
specular intensity, for larger momentum transfer, due to destructive interference within the 
rough layer. These are based on probability distributions for height, or slope, p(rj), etc., 
that require additional assumptions about surface properties. For instance, the commonly 
used Gaussian form = 6 _1 (27r)~ 1 / 2 exp[— £ 2 /(26 2 )], which is normalized and has second 
moment (£ 2 ) = b 2 , is not implied by the Gaussian form of f(S); only the second moment is 
common. While the correlations (4), (7) and higher are unambiguously determined by f(5) 
alone, all the higher moments for constitute additional assumptions. For real surfaces 
even b and w are usually not well known. We show below that from f(S) alone a Debye- 
Waller factor can be determined, but only up to order b 2 , if the perturbation approach 
for the rough wall is carried to second order. Other assumptions usually made to derive 
'Debye- Waller factors' will also be discussed. 

We will also point out that asymptotic aspects of roughness-induced scattering and ab- 
sorption can be obtained from mean values (b 2 , a 2 , etc.) only, independently of the details 
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of the models used (Gaussian or 1 Kq or similar). Going to second order perturbation is 
both necessary and sufficient for a full account of reflection and absorption up to quadratic 
terms (b 2 ). Higher perturbations are at least of order ft 4 for symmetrical roughness where 

/(*) = /("*)■ 



III. PERTURBATION APPROACH TO UCN INTERACTION WITH A ROUGH 
WALL 

As usual [4-6] a wall with micro-roughness may be divided into a volume V with an 
ideally smooth wall at z — 0, and a thin roughness volume V\ with partly positive and partly 
negative thickness £(p), where (£) = 0. For an incident plane wave ipi(r) = exp(ikj • r) = 
exp(— ikz cos 9i)exp(ikx sin 9i) approaching the surface from the vacuum side (z > 0) the 
neutron wave function ^(r) may be split into the unperturbed part ipo( r ) for the plane 
surface and a small perturbation ipi(r): 

V>(r) = Vo(r) + Vi(r) = Vo(r) + ^ 1} (r) + ^ 2) (r) (8) 
The first-order perturbation is obtained [4] as 

^(r) = -q J G(r\r')i> (r')d 3 r' (9) 

where the scattering potential q = Na(l — irj) = k^/Ajr = (A^ /47r)(l — irf) is determined 
by Na, a mean value for the number N of atoms times their bound-atom scattering length 
a. The imaginary part given by rj = —Im(q )/Re(q ) takes into account loss processes 
(nuclear capture, inelastic and, except for UCN, incoherent-elastic scattering), i] is small for 
wall materials of interest for UCN storage. k c0 denotes the critical wave number for total 
reflection at normal incidence on a flat wall of uniform scattering potential. 
The Green's function G(r\r') in (9) satisfies the equation 

V 2 G(r|r') + K 2 (r)G(r\r') = -4vr5(r - r') (10) 

where K(r) is the wave number: K = k in vacuum, and K = (k 2 — k^) 1 / 2 within the medium. 
For negligible refraction, G(r|r') = exo(ik\r — r'|)/|r — r'|. For the general expressions see 
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Appendix B. Asymptotic expressions for G(p, z\p', z') in the limit of large r have been used 



This method can readily be extended to obtain the second and higher order perturbations. 
To obtain the (n+ l) th order correction to ipi(r) we insert ^^(r) in the integral of (9), and 
below we will apply this method to extend roughness scattering to second order. Going to 
third or higher order would require correlations of higher than the second order. 

Instead of the 'distorted Green's function method' of [4] used here, an alternative Dis- 
torted Wave Born Approximation (DWBA) was used in [5]. In [6] it was shown that the two 
methods are equivalent, with one important caveat: A conventional 'Born approximation' 
uses the far- field Green's function expansion to obtain the factor exp[i(k — kj) • r'] inside the 
integral. This gives an adequate description only in first order. Second order calculations 
require also the near field of the Green's function, at z ~ 0, in the integral of (9) to obtain 
V>i 1} (r) at the rough surface. A variation of this method was used earlier in [7]. 

IV. GEOMETRY, BASICS AND MAIN RESULTS 

For a wave incident with a wave number k in the (zx) plane perpendicular to the wall 
at polar angle Q{ (measured from the wall normal), Fig. 1 shows projections of incident (i), 
reflected (r) and (elastically) scattered (s) wave vectors onto the (xy) plane. The scattered 
beam at solid angle (9, ip) is characterized by the in-plane momentum transfer parallel to 
the wall, q = (k s — k r ) par , where 



Since we consider only elastic scattering, the end point of q is restricted to the area inside 
the circle with radius k. For the part of ijj containing the outgoing wave and the evanescent 
wave inside the medium we use the expansion 



in [4]. 



q = k (sin9 — sin^) 2 + 4 sin # sin sin 2 



1/2 



(11) 



^(r)=^0r(r)+^(r) + ^(r) 



(12) 



where 




R(9i) exp(ikz cos 9i) exp(ikx sin 9i); z > 
S(9i) exp(Kiz) exp(ikx sin9i); z < 
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FIG. 1: For elastic wall interaction of UCN with wave number k we show the projections onto the 
surface plane of incident, (kj), mirror-reflected, (k r ), and diffusely scattered, (k s ), wave vectors. 
The in-plane wave vector transfer q is the component of k s - k r parallel to the surface. The vectors 
are bracketed to indicate projection onto the xy-pl&ne. For elastic scattering the end points of q 
are confined to the circle of radius k. 

is the unperturbed wave field for the flat wall, with 

, k COS 6; — IK; 

R ^ = I T~ 

K COS t>i + IKi 

and 

S(9i) = l + R(9i); m = yjkl - k 2 cos 2 6 t . 

^\r) and ip^ir) denote the first and second-order perturbation terms. 

Keeping terms ~ i] 1 in R, S, etc. we obtain the absorption corrections to the specular 
and scattered beams and, by comparison of net outgoing to incoming intensities, the loss 
coefficient for the rough wall. For details of the analysis see Appendix C. Here we summarize 
the main results for the intensities normalized to the incident flux. Taking the squared 
magnitude of (12) (outgoing waves only) results in cross terms labeled by (mn) where m,n 
= 0,1,2 stands for the order of perturbation: 

(a) (0,0): Specular reflectivity for the flat wall, derived from \^o\ 2 : 7( o) = \R{@i)\ 2 — 
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1 — /io(0i) where fio(Oi) = 2r]kcos(9i)/ K i0 is the absorption coefficient for incidence at 
angle 0* on a flat wall [7-9], with K i0 = (k% — k 2 cos 2 ^i) 1 / 2 . 

(b) (0,1): Interference term with first-order perturbation for specular reflection, arising 

from 2Re (VqV/^): 

V) = -^^b 2 (13) 

For r] = 0, there is no first-order interference. This shows that a 'Debye- Waller factor' 
for specular beam attenuation cannot be obtained from first-order perturbation. 

(c) (1,1): The first order intensity scattered into unit solid angle at (9, ip) is given by the 

outgoing flux with density |?/>^(r)| 2 , normalized to the incident flux, with the result 
(for UCN) 

(14) 



7 (11) = 4F cos Bi cos 2 6F{q) 
where 



, kcos9i kcosQ 

1-77 ^ + 



m = F( q ) = J A m^ s d 2 s (15) 

is the Fourier transform of the height-height correlation function (1) and represents 
the roughness spectrum [4]. For rj = 0, and for UCN (rather than cold neutrons in 
general) Eq. (14) agrees with Eq. (20) of [4]. Removing common constants, we write 
F(q) = [b 2 u?/(2ir)}L(q), where 

L(q) = exp(— q 2 w 2 /2) for the Gaussian model, and (16a) 

Of Ki \t(\ + 2o 2 7/; 2 W 2 l 

= n o 2 2M/2 ZL for the '^o-model' with t = 0.7709. (16b) 

(1 + 2q 2 w 2 ) 1 ' 2 K (t) 

From (14) we obtain the total diffusely scattered intensity, up to terms ~ b 2 and rj 1 : 
Pd = PD0 -2k 4 cos9 l [ d{lcos 2 9[fi (9, i )+fi (9)}F(q) (17) 

J(2tt) 



where 

Pdo = 4/c 4 cos#i 

'(2tt) 



Pm = 4k 4 cos 9i [ dtt cos 2 9F(q) (18) 

J(2n) 



is the total scattered intensity for rj = 0, and /jLq(9) = 2r]kcos9 / {k 2 Q — k 2 cos 2 9) l l 2 is 
the flat- wall loss coefficient for angle 9. For a « 1 the scattered intensity /(n) forms 
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a halo around the specular beam with small width of order a, as expected from ray 
optics. 

(d) (0,2): Interference of specular reflection with second-order perturbation, arising from 

2Re (V Vrf 2) ): 

I m = -Pdo[1 - A*o(0i)] + k 2 k 2 c0 cos9 l [ dQ^(9)F(q) (19) 

J(2tt) 

For r\ — 0, this term is required to satisfy unitarity. Up to order b 2 , the total outgoing 
intensity (specular plus diffuse) equals the incoming intensity. In other words, the 
'Debye- Waller' attenuation factor DWF for the specular beam equals 1 — p D0 , as it 
should. Carrying out the integration for p m we find that the DWF is close to but not 
identical to the factor exp(— 4b 2 k 2 cos 2 0,) ~ l—4b 2 k 2 cos 2 0j given in [5]. At this point 
it should also be mentioned (as is in [5]) that to arrive at the common Gaussian form 
for a DWF one has to make the drastic assumption that the wave inside the medium 
is given by the same function as the wave outside the medium, throughout the rough 
layer. This is not plausible for larger roughness, especially not for UCN, since only ip 
and its first derivative dip/dz are continuous at the surface. Using (13) and (19) we 
obtain the Debye- Waller factor for an absorbing, rough wall: DWF = |-R| 2 + /(oi) +-f(02) 
(up to terms ~ b 2 and ~ 77). 

(e) Combining the terms (0,0), (0,1), (1,1) and (0,2) we obtain the absorption coefficient 

for a rough wall: 

fj,(9i) = (incoming flux - outgoing flux)/(incoming flux) = 1 - (/(oo) + ^(01) + ^(02)) — 
= fx (0i)(l+k 2 c0 b 2 ) -2k 4 cos 9i [ dn cos 2 9[fi (9i) -fi (9)]F(q)-k 2 k 2 c0 cos 9i [ dVt^{9)F{q) 

J(2tt) J(2tt) 

(20) 

This agrees with the result derived in [7] using a somewhat different approach. Nu- 
merical results were not given in [7]. Higher-order perturbation would yield only terms 
~ b n with n > 2 . 
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V. APPROXIMATIONS AND NUMERICAL RESULTS 



Numerical evaluation of the double integrals over 9 and ip in (17)-(20) is greatly facilitated 
if one of the integrations can be performed analytically. For the Gaussian model, the ip- 
integration leads to the Bessel function I , but no analytical solution is known for the 
'ii'o-model'. Using the transformation k 2 dilcos9 = d 2 q = k 2 vdvd^, with q = ku, the ip- 
integration can be performed analytically for any roughness model since F(q) only depends 
on v, not on ip. It is evident from Fig. 1 that for v < 1 — s« (with Sj = sin#j), the 
•^-integration runs from — n to +7r. For 1 — Sj < v < 1 + Sj, the limits are ±ip u with 
sin 2 (-?/> u /2)) = [1 — (y — Si) 2 ]/(AvSi). Expressing cos# in terms of v and ip as cos 2 9 = 
1 — (si — v) 2 — Avsi sin 2 (-^/2), the -^-integrations in (17)-(20) can be performed analytically 
in terms of elliptic integrals, as shown in Appendix D. 

Before proceeding with numerical results, we point out the special case of fairly smooth 
surfaces with small mean-square slope (a 2 << 1) and a UCN wavelength small on the 
scale of the lateral correlation length w, viz. kw » 1. Smooth surfaces are expected to 
form when a special low-temperature oil 'LTF' [10] is sputtered onto a cold surface, then 
thermally cycled by slow liquefaction and re- freezing, as in [3]. Under these circumstances, 
F{y) is very small for v » (kw)^ 1 and the //-integration can be extended to oo. In this 
case, model-independent results are obtained as follows: The analytical expressions for the 
-0-integrals are expanded for small v as 1^ = a + a 2 v 2 + a 4 z/ 4 + ... and we use the identities 

poo 

271 / qdqF(q) = /(0) = b 2 (21) 
Jo 

POO 

2tt / q 3 dqF(q) = g(0) = a 2 (22) 
Jo 

POO 

2tt / q 5 dqF(q) = 4/i(0) = An 2 (23) 
Jo 

which follow from the properties of Fourier transform (15) and relations (4) and (7) between 
f(5),g(5) and h(5). 

In this way asymptotic expressions are obtained for p D0 and A/i(9i) / )io(9i) = \jJ,(0i) — 
W>(0i)]M>(0i): 

Pm (9 t ) = Ak 2 b 2 {c 2 - ^#7^ + 0[(M" 4 ]} (24) 
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and 

A/^VM^) = - 2(1 5 C2g2) {l " + ^(M" 2 ]} (25) 

where ( = k/k c0 < 1 and q = cos 6^. For perturbation theory to be valid, Pooi^i) mus t be 
<< 1 over the entire range of 0j. 

Taking only the first term in the curly bracket (p D0 46 2 A; 2 cos 2 Eq. (24) is con- 
sistent with the usual DWF with expansion (1 - Ab 2 k 2 cos 2 9i + ...), which is the standard 
Born approximation result of [5] (see also [12]). Note that for rj = the loss of specular 
intensity, 1-DWF, equals PD O (0i) [6], thus particle number is conserved. Figure 2 compares 
approximation (24) (first term only) with the exact numerical result for k c0 = 0.0802 nm _1 , 
a = 10~ 3 and w = 2 /im (b = aw/y/2 = 1.4 nm). The critical value k c0 is for the LTF oil in 
its glass state at ~ 110 K used in [3]. The agreement is excellent, indicating that for these 
parameters the Born approximation for the DWF is adequate. The numerical values of p D0 
for the Gaussian model are very similar to those of the .fTo-model used for Fig. 2. 

On a general basis, for a 2 « 1 the approximations (24-25) are valid over a wide range 
of 9i, but not near grazing incidence (6>j — > n/2) since in this range (1 — Si)kw is not » 1. 
Therefore the z/-range for integration over the full -^-range from — it to +7r (see Fig. 1) 
cannot be extended to oo. 

From (25) we deduce a, perhaps, unexpected result: over a significant range of 9i, namely 
where the curly bracket is positive (for ( << 1: 9i < 45°), the loss coefficient for the 
rough surface, is somewhat smaller than that for the flat wall, Ho(0i), with a relative 

difference of order a 2 . This contradicts a general statement in [7], p. 185, while a decrease 
was predicted for a certain model in [11]. The decrease found here for a « 1 is independent 
of the model and is plausible for a geometric-optics picture, where geometric reflections for 
the local surface orientation are incoherently superposed: For near-normal incidence the 
mean angle for incidence on inclined patches of a rough surface is always larger than for the 
flat surface, for which it is zero for Q, L = 0, and therefore the factor (A; 2 — A; 2 cos 2 (^)) ^ 2 
in Ho{6i) is reduced, on average. No reduction of Ho(9i) is seen for 'jagged' roughness with 
a of order 1. For the same values of a, etc. as for Fig. 2, Fig. 3 shows a comparison of 
numerical results for A//(#j)///o(#i) with approximation (25) for the K model. Again, the 
approximation clearly fails for grazing incidence (both (24) and (25) diverge as 6>j — > n/2) 
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FIG. 2: (color online) The total diffusely scattered intensity poo (from (18) for the Kq model) is 
compared to approximation (24) which is valid for fairly smooth roughness (a << 1). For the 
parameters used the agreement is excellent. 

but is an excellent representation at steeper incidence. The Gauss and K Q model give very 
similar results. 

For UCN confined in a trap the angular distribution is, in most cases, assumed to be close 
to isotropic (although deviations are critical, see below). To obtain directional averages of 
p D0 and A/i the region of grazing incidence that is not described by (24-25) has to be 
included, and it makes a significant contribution. The angular dependence of A/i(0j)/// o (0i) 
for the Gauss-model (not shown in Figs. 2 and 3) with the same parameters is close to that 
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FIG. 3: (color online) Angle dependence of roughness correction to the wall loss coefficient, 
Afi(9i) / no{0i) . The numerical result from (20) for the Kq model is compared to approximation 
(25). The agreement is very good except for 6i near 90° (glancing incidence). In this plot the scale 
changes from n = 3 (expanded scale) for 0, < 70° to n = 1 (direct scale) for 6i > 70°. 

for the i^o-model and the isotropic mean values 

(Z) = 2 I c.dciZie,) (26) 
Jo 

agree within 15 percent. Fig. 4 shows the ^-dependence of (A/x) / (/j,q) for both models. The 
reference average for the flat wall is given as [7-9] 

4r/C/3 for ( « 1 
ITT] for C = 1 



( Ai0 ) = 2 ?? [ a rcsmC-C(l-C 2 ) 1/2 ]/C 2 ^ 
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FIG. 4: (color online) Mean value of roughness correction, (A/z(0j)), for isotropic UCN flux, referred 
to the flat-wall average (/io(#i)) given in (27). For the parameters given, the numerical result for 
both models are compared to the semi-analytical approximation (28). Both for the Gauss and the 
Kq model, the agreement is reasonably good, with maximal deviations of ~ 20% over a wide range 
of parameters. 

The main features of (A//) / (fi ) are a sharp peak of height (bk c o) 2 at ( = with half- 
width A£ ~ \/2/(wk c0 ) and a drop-off approximately ~ C~ 3 ^ 2 - A reasonable fit, with 
maximum deviations of ~ 20% over a wide range of k c0 and of (small) a and b and (large) 
w, both to the Gauss- and K$ micro-roughness models, is 

<M - A; c > 2 [1 + 3(C W2) 2 ] - 3/4 (l - ^C 4 ) (28) 
which is shown as the dashed curve in Fig. 4. The isotropic-flux average of diffuse fraction 
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Pdo is 

(pm) = 2k 2 b 2 (29) 

where we used the one-term approximation of (24). 

In applications to actual UCN storage, approximations (28,29) can be used to determine 
further averages analytically. For a monochromatic UCN spectrum under gravity the aver- 
aging is over wall-interaction height h from the trap bottom to the 'roof or to the maximum 
height reachable for a given energy (ho, in units of maximum jump height). Further aver- 
aging over the spectrum (over ho) is usually also required. For instance, the sharp peak of 
(A/x(0j)) / (fio(0i)) at C = implies a large loss rate for UCN in contact with the wall near 
their maximum jump height h where ( — > 0. Taking into account that the flux [~ (h — h)dS\ 
incident on surface element dS vanishes as h — > ho, a detailed analysis is required. It shows 
that the loss enhancement strongly depends on the trap geometry and UCN energy and is 
significant for very low energy UCN 'hopping' in small jumps on an essentially flat horizontal 
surface, for instance the rim of a cylinder with horizontal axis, as in [3]. 

VI. DETAILED BALANCE REQUIREMENT AND SIMPLIFICATIONS 

A trapped UCN gas will acquire or maintain an equilibrium isotropic distribution only 
if its interaction with a rough wall (or magnetic field irregularities) satisfies the detailed 
balance requirement discussed in [7], p. 96. The flux undergoing scattering from solid angle 
Qj = ifi) to Q = (9, ip) must equal the flux for the reverse process Q — > Qj, in accordance 
with the fundamental principles of time-reversal invariance and micro-reversibility. Since the 
flux incident at Qi (or Q) is proportional to cos 0j (respectively cos 9) any diffuse scattering 
distribution I sc (fli — > Q) must satisfy the symmetry requirement 

i sc (Qi ->■ n) = ifa, n) cos e (30) 

where fl) is a function symmetric in (9i, 9) and (ipi, tp) and cos 6? is the Lambert factor. 
The scattering distribution Iny of Eq. (14) satisfies this requirement. 

Various simplifications of the micro-roughness scattering distribution (14) (for r] — 0) have 
been used in [13-15]. The limit of maximal diffusivity is reached for a "dense roughness" 
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model where kw — > 0, hence a ~ b/w — > oo. This leads to [14,15] 

I SC ^B [cos 6i cos 9} cos (31) 

with _B = 2k 4 b 2 w 2 /ir and a diffuse fraction p D0 = (2ir/3)B cos ^ that depends on incident 
angle #j and UCN energy (~ k 2 ). It has been pointed out in [15] that for kb at most of order 
1, as required for perturbation theory to be valid [4], the limit kw — > 0, hence B — > 0, is 
strictly justifiable only for very small diffuse fraction p D0 . 

As a further simplification, averaging over incident energy may be performed but could 
be a coarse approximation for broad spectra, as for confined UCN. To arrive at the simplest 
form of (30), with I = const., requires further averaging of [cos 9i cos 9] in (31). Assuming 
isotropic distributions in 9i and 9 we obtain (cos 9i cos 9) = 4/9. Thus, (31) becomes 

(J sc ) = (l)cos9 (32) 

with (i^ = 8 (A; 4 ) b 2 w 2 /{9n) and a total diffuse fraction (p D0 ) =*(l)- Like (31), Eq. (32) 
is justified only for small (pdo)- 

In the "dense roughness" limit the loss coefficient fx(9i) of (20) can be averaged analytically 
[7] with the result (//) = (fx ) (1 + k 2 Q b 2 ), where (/x ) is given in (27). However, except for 
very small (pno) the correction term k 2 b 2 = k* Q (p m ) / ((k A ) w 2 ) would become large as 
k c oW — > 0, conflicting with the perturbation theory requirement. 

It should also be noted that in averaging over 9i we lose the proportionality of I sc to cos 
which implies a small scattering probability for glancing incidence (9 { — > ir/2). For UCN 
stored in traps with high geometrical symmetry this dependence can give rise to almost sta- 
tionary orbits for UCN "sliding" along a concave surface (see Sec. 8). For other geometries, 
only a few reflections, or where these details are not important, as in the computer code 
tests of [13], (32) is a useful, simple approximation satisfying detailed balance. 

As an example of a roughness model that does not satisfy detailed balance we mention the 
scattering distribution (E6) derived in Appendix E for the macroscopic, ray optics limit. In 
this symmetrical in (#;, 9) and (ipi, ip), lacking the extra factor cos^. This violation 

of detailed balance can be tolerated if only a few reflections have to be considered, as in a 
neutron guide [4], but attempting to use this model for thousands of consecutive reflections 



16 



in a Monte Carlo simulation, as required for long UCN storage (see below), results in the 
loss of isotropic equilibrium. We observed that UCN accumulated in grazing-angle orbits, 
depleting other regions of phase space. 

VII. MONTE CARLO SIMULATION TECHNIQUES 

In Monte Carlo simulations of UCN propagation and storage, as in [7,14,15], two comple- 
mentary ways of implementing a given scattering probability distribution P(9,(p), e.g. Eq. 
(14), have been discussed and used. 

One ([7], p. 98) is based on mapping the joint probability distribution P(9,ip) onto the 
ranges to 1 of uniformly distributed random numbers (r.n.) x, y, z, etc. For instance, 
for the 'dense roughness' limit of Eq. (14) used in [14] in the form P(9,Lp) ~ cos 9i cos 2 9, 
we can determine 9 and <p from x, y, setting cos 3 # = x and <p/(2ir) — y — \- In general, 
the mapping procedure involves indefinite integrals of P(9, 0) and becomes very tedious if 
numerical integrations are required as for (14). In the simple case P ~ cos 2 9 the indefinite 
integral J cos 2 9d(cos 9) is ~ cos 3 9, hence the prescription cos 3 9 = x. Rendering this method 
efficient requires 3D tabulation where the tables need only be interpolated at each UCN 
reflection during a simulation run. As a further simplification, once 9 has been determined, 
if is found faster from the conditional probability of ip, given 9: p(tp) = P{9, ip) / Pg(9), where 
Pe{9) is the integral probability for 9. For the Gaussian model the numerical work involves 
only the incomplete Bessel function, which is readily available in computer code. 

For our simulation (below) we chose the second, often faster method used in [15], to 
determine 9 and in a way consistent with the probability distribution. First read from a 
table the total diffuse-scattering probability p m for given incidence parameters k and 9i. A 
r.n. xq then determines whether the reflection is specular (if x > p m ) or diffuse (x < p D0 ). 
If diffuse, a trial set (9, ip) is obtained as cos# = x and (p = 2n(y — \). Next, find the 
probability P(9,(p) for the trial set and compare it to a fourth r.n., z. The angles (9,p) 
are accepted if z < P{9, ip)/ P max , where P m ax is the maximum of P(9,ip) for incidence at 
k and If z > P(9,<p)/P max , the process is repeated with new r.n.'s x,y,z. This scheme 
can be greatly abbreviated if 1-nP max < 1. In this case rolling three dice (x, y, z) only once 
is sufficient, x and y determine 9 and <p, as before. If z > p D0 , the reflection is specular, 
otherwise diffuse at angles (9, <p). It can be verified that this faster scheme also respects the 
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probability distributions correctly (as long as 2nP max < 1). 

In all these schemes, wall losses and beta-decay are irrelevant for the choice of reflection 
angles. However, we keep track of the accumulation of net loss at each reflection and use 
the overall loss factor as a weighting factor for UCN that have survived losses and can be 
counted in a detector. The loss factor contained the second-order correction of Eq. (20) and 
was implemented in tabulated form. 

VIII. MONTE CARLO SIMULATION FOR NEUTRON LIFETIME EXPERI- 
MENT [3] 

The neutron lifetime experiment [3] reports a value r n = 878.5 ± 0.7 sta t ± 0.3 SJ/S s which 
is 7.2 s (or 0.8%, or > 6 standard deviations) away from the current particle data world 
average 885.7 ±0.8 s [2]. The quoted precision is better than for any previous single lifetime 
measurement. Since data interpretation for this experiment relied heavily on computer 
simulations we have performed independent simulations using the roughness model outlined 
above. In this system [16,3], shown in Fig. 5, a cylindrical or a 'quasi-spherical' vessel with 
an opening rotates about a horizontal axis during the various steps of a cycle: Filling with 
UCN with the hole pointing straight down; rotation to a 'monitoring position' for spectral 
cleaning; storage for a 'short' time (300 s) or a 'long time' (2000 s) with the hole in the 
upright position; then emptying in 5 steps at intermediary angles. This scheme provides for 
a spectral analysis, with a batch of higher-energy UCN counted first and the slowest UCN 
pouring out at the last stage while the opening moves to the vertical down position. The 
detailed measuring scheme [3] is complex, involving not only energy dependent data for the 
5 energy bins but also traps with different mean free path for the UCN. This combination 
allows, both 'energy extrapolation' and 'size extrapolation' to the neutron lifetime. 

Data analysis is also complex, relying heavily on computer simulation to determine the 
parameter (7) = (nf) jr\ against which measured inverse storage lifetimes are plotted for the 
extrapolation. The average (fif), with wall collision frequency /, is the wall loss rate (mean 
loss per second). 

For a narrow cylindrical trap the authors of [3] compared simulations with the measured 
time spectra for the first counting interval following short storage (300s) (see Fig. 14, lower 
part, of [3]). They concluded that a certain minimum roughness (diffuse fraction > 1%) of 
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FIG. 5: Schematic of the apparatus used in [3]. UCN from a guide tube 1 pass through a large 
entrance/exit channel to a cylindrical or 'quasi-spherical' storage chamber 8 that can be rotated 
about a horizontal axis. After a short or long holding time, UCN in five spectral ranges are, 
successively, discharged through the same channel to the detector 12. The system 14 is used to 
coat the walls with a special low-temperature oil. 

the trap wall was required to obtain agreement. This is a crucial step with respect to the 
neutron lifetime, and therefore we have performed independent simulations of short and long 
storage cycles for a narrow cylindrical vessel with radius 38cm, width 14cm, and aperture 
62° for the opening, to investigate this issue using the near-mirror reflection model with 
long-range roughness (a « 1) described above. 

In our simulation UCN are generated at the beginning of the monitoring phase (at trap 
angle 30°) over a horizontal area inside the trap volume and near its bottom. The distribution 
is isotropic with polar angle 9, referred to the up and down directions, derived from a r.n. 
x as cos 6 = x 1 / 2 . The energy distribution is the same as used in [3], namely the Maxwell 
distribution modified by a heuristic attenuation factor exp(—ho/h sp ) with h sp = 0.6 m (not 
specified in [3]). We used a maximum (minimum) energy of 0.95 m (0.05 m). 

The time sequence of monitoring, storage and emptying in five steps was the same as de- 
scribed in [3] . The wall reflections were assumed elastic except for the (very small) Doppler 
shift for wall reflections during times when the trap rotates at 5.1°/s (respectively 17°/s 
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on the way from 0° to 40° for a rotation time of 2.3s [3]) and the reflection is diffuse. For 
perfectly cylindrical geometry the overall shift is very small. Reflection points were calcu- 
lated analytically from the previous reflection position, energy and take-off angle, assuming 
perfectly cylindrical geometry. At each reflection, the wall loss was incremented according 
to the loss coefficient /i for this reflection. A run was terminated at the time when the 
UCN passed through the opening without immediately falling back into the trap. For the 
purpose of this work the exit channel could be neglected and those UCN that had escaped 
loss, as determined by the integrated loss including beta-decay, were counted into time bins 
of 10s or Is resolution. For consistency with [3] we used the quoted values rj = 2 x 10~ 6 
and r n = 878.5 s in the simulation. The bins were combined to the five counting intervals, 
and short runs with storage time 300 s were compared to long runs (2000 s) as in the actual 
experiment to obtain the storage lifetime as a function of mean values (7) that were also 
calculated for each counting interval. 

Figure 6 shows simulated count rates, in 20 s bins, versus cycle time t, starting from 
'monitoring' for 300 s where the trap is at position 30° tilted away from vertical. The data 
are normalized to 10 6 runs (one run for each UCN started at t — 0). Data are shown for 
both short and long storage and for both the Gauss and K models, for a fairly smooth 
surface with similar parameters, as shown in the inset. They correspond to a — 7.1 x 10~ 4 
for the Gauss-model, and a = 2.4 x 10" 4 for K . An example of steeper roughness with 
a = 0.7 (for b = 2nm and w = 4nm) is shown for comparison. 

We will focus on two features of Fig. 6, and their implications: first the time spectra for 
the first counting interval following short storage; then the count-rates during long storage. 

First counting peak 

To allow a comparison with the time spectra for the first peak measured in [3], we show 
in Fig. 7, at a higher resolution of 1 to 2s, an expanded view of the interval from 600 to 750s 
(on the time scale of Fig. 14 of [3] this corresponds to the interval 740 to 890s). The data 
of [3] are included in Fig. 7 as curves 1 and 2. Prior to discussion we point out differences 
between the two simulation approaches. 

The authors of Ref. [3] used a roughness model characterized by a single number, the 
diffuse fraction for the trap walls, but did not state which roughness model was used. The 
simulation included the exit channel (designated secondary volume and UCN guide in [3]) 
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FIG. 6: (color online) Simulated count rates, in 20 s bins, versus cycle time t, starting from 
'monitoring' for 300 s where the trap is at hold position 30° titled away from vertical. The data 
are normalized to 10 6 runs (one run for each UCN started at t = 0). Data are shown for both short 
and long storage and for both the Gauss and Kq models for a fairly smooth surface with similar 
parameters, as shown in the inset. They correspond to a = 7.1 x 10 -4 for the Gauss- model, and 
a = 2.4 x 10 -4 for Kq. An example of steeper roughness (a = 0.7) is shown for comparison. We 
plot the simulated counts N plus one to allow a distinction between N = and N = 1 on the log 
scale. In the simulations of Fig. 6 (only) we used a constant angular speed of 4.3° /s for all trap 
rotations. 

but the authors did not describe the geometrical model used for this complex structure. It 
included (see Fig. 5) the trap exterior and interior (which some UCN temporarily re-enter on 
their way to the detector), the wide conical channel section, the valve selector with its gaps, 
the curved guide and the detector with its window that accumulated an oil film as a result 
of repeated surface coating using the evaporator 14 shown in Fig. 5. The surface properties, 
including roughness, of the oil-coated channel walls with their temperature gradient, varying 
from the trap temperature to near-room temperature for the detector, are not discussed. 
Our simulation did not include the exit channel. We used the roughness models and 
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FIG. 7: (color online) Comparison of measured [3] and simulated time spectra for the first counting 
period following short storage. For clarity, part (b) shows an expanded view of the initial 60s of 
the full 150s interval shown in part(a). Some of the simulation curves of part (b) are also included 
in (a). For a discussion see the text. 
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parameters listed for curves 3-12 of Fig. 7. In Fig. 7b the initial part is expanded for clarity. 
For these counting spectra following short storage the details of wall loss are insignificant and 
the results using (20) are indistinguishable from those for the elementary, flat-wall expression 
that was used in [3]. (During the long holding period of 2000s the details of loss coefficient 
may become important, depending on the roughness parameters.) 

Curve 1 represents the measured data (Fig. 14 of [3]), using a scale approximately ad- 
justed to counts per 2s for 10 6 initial UCN. Curve 1 practically coincides with the Ref. [3] 
simulations for 10% and 100% diffusivity. Curve 2 is the simulation, in [3], for 0.1% diffusiv- 
ity Curve 3 is our simulation (for the trap only) for an ideally smooth surface. For curves 
4-9 we used the simplest scattering distribution consistent with detailed balance, Eq. (32), 
with diffuse fractions 0.03, 0.1, 0.3, 1, 10 and 100%. For consistency with [3], we included 
large values of (pdo), disregarding the restriction outlined following Eq. (31). Curves 10-12 
are based on the micro-roughness result (14) (for rj = 0) and the i^ -model. The parameters 
for curves 10-12 represent increasingly steeper roughness (a = 2.4 x 10" 4 , 0.71, and 1.4) and 
increasing (p D0 ) (0.16, 0.55, and 11.3%, where the last value is at the limit of perturbation 
theory validity). 

Fig. 7 shows the large spectral distortion and time delay due to the exit channel. The peak 
widths at half maximum for curves 1 and 2 are 3 to 6 times larger than for the simulations 
excluding the exit channel. Thus the channel response function to the relatively short pulse 
of UCN pouring out of the trap essentially determined the measured time spectra. 

For consistency with [3] we focus first on the simulations based on (32). The minimum 
half-width At=10s is observed for (p D0 )=0.1%. Both for (p D0 ) >1% and (p D0 ) « 0.1% 
(including zero), the curves are ~2 times wider, while the curves for 0.03% and 0.3% are 
intermediary. The difference between 0.1% and >1% diffusivity makes plausible the (much 
smaller) difference of broadened peaks between the Ref. [3] simulations (including the chan- 
nel) for (p D0 )=Q.l% (curve 2) and 10 or 100% (curve 1). Comparing to the measurement 
the authors rejected the possibility of diffusivity << 1% for the trap. However, curve 3, for 
{Pd Q )=0, has about the same larger width, namely Ai=23s, as for (p DQ ) >1% (curves 7-9). 

The same type of broad line is seen for curve 10 that is based on Eq. (18) for small 
a. The intermediate (steep) roughness model used for curve 11 (12) results in intermediary 
(broad-type) linewidths. 

The difficulty of distinguishing between roughness models for the trap on the basis of 
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measured curves broadened 3 to 6 times by the response, or resolution, function is exac- 
erbated by two aspects: First, by the uncertainty introduced by any simplification of the 
complex channel geometry, as needed to make the simulations feasible. Secondly, in this 
experiment (and perhaps any UCN experiment conducted so far) neither the primary source 
intensity and spectrum nor the transmission characteristics of connecting guide tubes, shut- 
ters and other components, here including the channel that is passed by the UCN also during 
trap loading, is known well in absolute terms. Thus, the area under the counting peak can 
practically not be used as a criterion for adoption or rejection of a roughness model for the 
trap. In the present case, the observed variation of peak area in Fig. 7 is a reflection of 
spectral mixing between the five spectral intervals, as discussed below. 

We conclude from this peak analysis that a "soft-roughness" model with a « 1, as is ex- 
pected for the temperature-cycled [3] oil-glass coating used, can, probably, not be excluded. 
This has the consequences outlined below. 

'Spill-over' during long storage 

Returning to Fig. 6, we draw special attention to the finite count rate that appears 
during the storage intervals, t = 300 to 600 s for short storage and t = 300 to 2300 s for 
long storage. At the end of monitoring at trap angle 30° most UCN with energies exceeding 
the height barrier (~ 56.5 cm above trap bottom) have spilled out. However, some higher 
energy UCN move, preferentially, along the cylinder axis and little in other directions. In 
this highly symmetric trap geometry they have survived in nearly stationary orbits since, 
for a « 1, the reflections are nearly mirror-like with small reorientations of order a if 
the reflection is diffuse. Thus the random walk toward isotropy is very slow. As the trap 
rotates to angle 0° at t = 300 s the barrier rises by ~ 14 cm to 70.5 cm from the bottom, 
and the count-rate in Fig. 6 drops to zero. However, since UCN with energies exceeding 
this new barrier still remain, the same slow relaxation toward isotropy gradually results in 
steeper orbits. At t ~ 600 s some UCN have reached the new barrier and we see a fairly 
constant 'background' count rate throughout the remainder of the long storage interval up 
to t = 2300 s. The range of parameters a, b and w showing this general tendency is fairly 
broad, and there is no qualitative difference between Gauss and K models. We will show 
that this 'background' could be overlooked in actual measurements. The same slow diffusion 
in phase space is also expected to mix the spectral regions for the five counting intervals, 
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and the net result of slow relaxation could give rise to a significant systematic error. 

For comparison with experiment [3] we note that the counts in Fig. 6, for 10 6 runs, 
correspond to ~ 12 double cycles (one cycle each for short and long storage). The instrument 
background of ~0.02 s" 1 (Fig. 10 of [3]) is equivalent to ~ 5 counts/20 s in the region 600 s < 
t < 2300 s of extra holding time for long storage in Fig. 6. This interval has no counts without 
the relaxation process since the instrument background was not included in our simulations. 
Thus, the UCN 'spilling over the threshold' during this time (~ 0.8/20 s, or 67 in total for 
the K model) are missing from the UCN count, Ni mig , following long storage, but are not 
missing from the short-storage counts N s h or t- In the actual experiment the relative change of 
measured 'background' (0.8/5 ~ 15%) is below the significant background variations shown 
in Fig. 10 of [3] and would hardly be noticed. However, calculating from N s hort/ Ni ong an 
inverse storage lifetime, A = l/r st = [ln(N short / N long )}/ At, with At = ti ong - t short = 2300 s 
-600 s = 1700 s, disregards the 'spill-over' and would result in a significant uncertainty of A, 
as shown in the following estimation. 

We assume that the spill-over rate n S pui(t) originates exclusively from the highest-energy 
interval #1 (of the 5 intervals) and that entry into (exit from) spectral interval #1 by angular 
diffusion from (into) #2, or any other interval, is negligible. Thus the number Ni(t) of UCN 
in #1 decreases during long storage as 

1 dN, 1 1 

-n sp iii{t) (33) 



JVi(f) dt r st N^t) 

In [3] r st is close to the lifetime r n , so we can approximate N\(t) by Ni(t s hort) exp[— (t — 
t short) An)] in the last term of (33). Integrating (33) from t short to t long and equating N^tshort) 
with the counts for interval #1, Nshort, and Ni(ti ong ) with Ni^ ong we obtain 

\ — _L ^ _L l n Nl,short _ C S piU 

r st At N lt i ong N hshort At 

The first term on the right is the A-value disregarding spill-over. The second term is the 
correction with 

C spui = ^2,n spi ii(ti) exp[(tj - t short ) /r n ] (35) 

i 

where the sum is over the count rates n sp ui(ti) from t s h ort to ti ong , with a weight factor 
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exp[(ij — t s hort) /T n ]. Finally, the correction becomes 
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(36) 



This value is given in the inset of Fig. 6 and amounts to significant corrections to the 
inverse storage lifetime for interval #1. If inter- interval diffusion is included the values for 
other spectral intervals, and presumably the energy extrapolation as well, are affected to 
a significant extent. On the other hand, the measuring scheme of [3] is more complex, 
involving also the larger traps including a 'quasi-spherical' vessel, which in the simulations 
was replaced by a cylinder in a way not discussed in [3]. We will point out below that our 
simulation for a fairly smooth surface of a narrow cylindrical trap leads to questions also 
regarding the ^-values of the extrapolation, namely the mean values (7) for the spectral 
intervals 1-5. 

For the example with steeper roughness in Fig. 6, relaxation is faster and the large 'back- 
ground' in the region 300 s < t < 600 s might be visible in the experimental data, in spite 
of blurring by the exit channel. 

Mean 7- values 

For the narrow cylindrical vessel, Fig. 8 combines (7) values for the five spectral intervals: 
on the x-axis for an 'intermediate' roughness with b = 2 nm and w = 4 nm, and on the 
y-axis for a 'smooth' roughness with b = 0.5 nm and w = 3 /iin. The data for short and 
long storage are plotted separately and the differences are seen to be substantial. They are 
not explained by the small spectral cooling during storage due to larger wall losses for higher 
UCN energies. Spectral mixing and spill-over are more important. The large differences are 
also seen in the mean collision frequencies (not shown) for the five intervals. 

Mean spectral energies 

Spectral mixing is also seen in the 'erratic behavior' of mean spectral energies (height (h) 
and difference, A (h), between short and long storage). In Fig. 9, A (h) / ((h)) is plotted 
versus ((7)) (averaged over long- and short storage) for the 'smooth' roughness parameters. 
A spectral change between short and long storage implies a change of detection efficiency 
that is not discussed in [3]. While the difference in mean UCN velocity is reduced by 
the energy boost through the fall height of ~ 1.5 m between storage vessel and detector, 
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FIG. 8: (color online) Simulation results for the mean value (7) = (fj) /rj, for the five spectral 
ranges 1 to 5 of [3], plotted for short storage (300 s, solid lines) and for long storage (2000 s, 
dashed lines). For the Kq model, the abscissa shows an example of 'intermediary roughness' with 
6 = 2 nm, w = 4 nm, thus a = 0.7. The ordinate values are for b = 0.5 nm, w = 3 nm, thus 
a = 2.4 x 10 -4 . The seemingly 'erratic' behavior in either case may be explained by mixing between 
the five spectral intervals and loss due to 'spill-over' during the storage time as discussed in section 
8. The statistical errors for the (7) values are insignificantly small. 

the transparency of the oil coated detector window for an isotropic UCN beam is a steep 
function of velocity in the region of interest, with roughly a change of 4% for 1% of UCN 
energy change. As a result the UCN spend more time in the entrance/exit channel and have 
a greater chance of being lost, for instance through gaps. 
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FIG. 9: The 'cooling' effect during long storage versus short storage. Instead of a small systematic 
drop (Ah) of mean spectral height (h) as a function of time, expected since the wall loss increases 
with UCN energy, the simulation shows fairly large and 'erratic' values of (Ah) for the five spectral 
ranges. Although there is cooling for the full spectrum, the lowest-energy intervals show negative 
values of (Ah) ('heating'). The simulation data are for the Kq model with fairly 'soft roughness' 
(a = 2.4 x 10~ 4 ) and are explained by inter-spectral mixing and by 'spill-over' (see section 8). The 
((7)) and ((h)) values are averages over the spectral ranges and over the two storage times (300 s 
and 2000 s). 

IX. CONCLUSIONS 

Using models of microscopic roughness characterized by height-height correlation func- 
tions of the Gauss and j Kq type (for solid vs. liquid surfaces) we have carried perturbation 
calculations to second order, recovering a previous result [7] for roughness-induced wall loss 
for UCN that has, apparently, found no attention so far. It may be surprising that under 
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certain conditions roughness reduces the wall loss probability. Furthermore, the second- 
order approach provides absorption corrections to the diffuse scattering distribution as well 
as the leading terms of a 'Debye- Waller factor for roughness' that is consistent with previous 
derivations but does not require the additional assumption of a height probability distribu- 
tion function made, e.g., in [1] and [5]. It is shown that for a fairly smooth wall with small 
mean-square slope the results are model independent and are determined only by the mean 
values, mainly for height (b 2 ) and slope (a 2 ). 

As the central application of roughness scattering we describe independent simulations 
relating to neutron lifetime experiment [3] which reported a significant deviation from the 
world average and very high precision. The possibility of a fairly smooth surface for the 
temperature cycled oil used as the wall coating in this experiment, in connection with the 
possibility of nearly stationary orbits in any highly symmetric trap geometry, raises questions 
about the reliability of the extrapolation method used to extract the lifetime value from the 
storage data. 
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Appendix A: Correlation functions 



Application of the Laplace operator twice in Eq. (7) yields for the Gauss model: 



h(5) = (K(pMp + 5)) = 4 ( 1 - ^ + ^ ) eMS 2 /(2w 2 )} 



(Al) 



with mean square curvature 



2 « 2 2b " 



Using properties of the modified Bessel functions K n (with n = 0, 1,2...) we obtain for the 
K model. 

h(5) = k\ 



2 Q(z) 

hK Q(t) 



(A2) 



with 



= [(5 2 + 5 2 )/(2w 2 )] 1/2 ; t = 5 /(wV2) 



Q{z) = 8 



Ko(z) 



5 2 + 5 2 



.K«(z) 



and 



4 = Mo) = f£ 



= 2.103a 2 /w 2 



(A3) 



(A4) 



[Ki(t) 

The expressions for the mean squared curvature hold for the choice t = 0.7709, for which 
the two models have identical mean-square slope a 2 . This value of t is the solution of 
K 2 (t) = 5K (t). 
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Appendix B: Green's function for a flat surface 



Eq. (10) describes the wave, with wave number k, generated at point r = (p, z) by a 
point source at r' = (p',z') near the surface z = of a semi-infinite slab of material with 
scattering potential k 2 /4iT. The solution may be represented by a Fourier expansion in terms 
of in-plane waves exp(ik|| ■ p) as [17] 

G(p, z\p', Z ') = -^J d\^»-^g{z\z') (Bl) 

The one-dimensional Green's function for k z satisfies the equation 



d 2 g(z\z>) 
dz 2 



+ K 2 (z)g(z\z') = -4n5(z-z') (B2) 



where the wave number component perpendicular to the wall is k z = {k 2 — k 2 ) 1 / 2 , and 
K z = k z in vacuum and K z = k' z = (k 2 — k 2 ) 1 ^ 2 inside the medium. For UCN with k 2 < k 2 , 
k' z = in = i(k 2 — k 2 ) 1 / 2 is imaginary apart from a very small real part due to the loss 
contribution %r\ in k 2 = k 2 (l — it]). The solution of (B2) for outgoing plane waves can be 
written, for z > z', z > 0: 

g(z\z') = ^ . e ikzZ e KZ ' 

4-7T7 7 

" —e lk * z (l + Kz') (B3) 



and for z < z', z < 0: 

g( z \ z ') = — e ^[e^ z '-Re KZ ') 

K 

—e KZ {l+tk z z') (B4) 

with R — (k z — in) / (k z + in) . 

(B3) and (B4) were derived assuming z' < 0. The two-term approximations in (B3) and 
(B4) are valid for 'micro-roughness', where k c0 b « 1, and they hold for positive or negative 
values of z' within the rough layer. This is a direct consequence of the continuity of ip and 
dip /dz at the interface. We note that the first two expansion terms are sufficient for a full 
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description of UCN interaction with a rough wall up to terms quadratic in the roughness 
amplitude b. 

The same kind of approximation can be made for the plane-surface wave ipo(r') in the 
integral of Eq. (9). Thus, using (B3) the z'- integral in (9) from z' = up to the roughness 
amplitude across the rough layer is proportional to £ + (k + /tj)£ 2 /2 where Ki = 

{k 2 c -k 2 cos 2 Qifl 2 . 

The full Fourier expansion (Bl) is needed in (9) in second-order perturbation where the 
first-order perturbation ■0i^( r ) obtained for z « in (9) is subsequently inserted in the 
integral with if) replaced by to obtain i[)^\r). 

For flux calculations, we need the far field expressions for G(p, z\p' , z') for z > 0. Using 
the two-term approximation for its z'-dependence (the last bracket in (B3)), Eq. (15) of ref. 
[4] reads 

9 h p ikr 

GW > z,) -v^— e ~* n (B5) 
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Appendix C: perturbation calculation up to second order 

Derivation of the results listed in section 4 is lengthy and we will only summarize a 
crucial steps, using the same symbols as in the text. 

(a) Integrations over p: 

By definition of £(p): 





By definition (1): 




By definition (15): 




£(p)£(p'y q - (p - p Vprfy = A \ f(5)e^ s d 2 5 



(27r) 2 AF(q) = (2ir) 2 AF(q) 



(b) UCNflux: 



Incoming flux for velocity vucn = (with m = neutron mass): 



in 



Av UC N cos 9i 



Outgoing flux into solid angle dVt: 



$out = v UC N\^out\ 2 r 2 dfl 



For the specular beam at angle Q r (9 — ip — 0) 



ikz cos 6i ikx sin 0i 
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we can use the plane wave expansion in spherical harmonics [18] for large distance r: 



o ikr °° + ; o p ikr 

^ r -r £ £ )^(^ *0 = T — ^ - Q r) (C6) 

2=0 m=-J 

This allows simple integration of the interference integrals. 
Loss terms: 

On each step the terms containing the loss coefficient rj are developed to first order in 
rj. For instance, for k z not too close to k c0 : 

«=V^ r *2 = «o-?^ (C7) 

2 /to 

with k 2 c = k 2 (l - ir?) and k = a/ k 2 - k 2 z 
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Appendix D: Integrations in q-space 



Using the transformation dVt cos 9 = vdvdip and the relationship 

cos 2 9 = 1 - (si - v) 2 - Avsi sin 2 (VV2), 

the integrations over ip in (17)-(20) can be performed analytically. The integrations run 
from to ip u where 

• J 1; f° r v<l-s 

sm — =< i — — (Dl) 

Since the -^-integrals needed in (18) to (20) do not appear to be readily available from tables 
or websites we list two results 

^ | 2v 1/2 E(m 1 ); v < 1 - s t 

coaOdip = I V ' (D2) 

y 2a 1 l 2 E(m^ 1 ) + 2v + a- 1 / 2 K(m^ 1 ); 1 + Si >v>l - Si 

and 

^ 2C 2 cos 2 9 - 1 , ( 1 2 M ; 1/2 K(m 2 ) - Au^E^); v < 1 

— dip — < 

o V 1 -C 2 cos 2 ^ y 2u~ l ' 2 F{u l \m 2 ) - AufE^m^) + 4C 2 (-^+) 1/2 ; 1 + s< > i/ > 1 

(D3) 

where [19] 

F(u\m)= / (1 -msin 2 /3)~ 1/2 rf/3 (D4) 
Jo 

is the incomplete elliptic integral of the first kind, 

E(u\m)= (1 -m sin 2 /3) 1/2 d/3 (D5) 
Jo 

is the incomplete integral of the second kind, and K(m) = F(||m) and E(m) = E(||m) are 
the corresponding complete elliptic integrals of the first and second kind, respectively. The 
arguments in (D2) and (D3) are 

C = v ± = 1 - (y ± Sl ) 2 ; «± = 1 - CV; a = 4z/ Sj ; m, = ±; m 2 = sin Wl = (^±)V2 
Compared to numerical double integration, the use of these readily available functions in 
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one of the integrals as well as tabulation of the numerical results for use in simulations is of 
enormous benefit. 
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Appendix E: Macroscopic limit 



In the geometrical optics model, reflection is considered as an incoherent superposition of 
rays reflected at each speck of a fairly smooth surface (a « 1) as if it were a plane inclined 
at the local slope. The analysis is based on a slope distribution function p(x) that has to be 
postulated independently of the correlation function g(S). For instance, p(x) does not have 
to be Gaussian if g(8) has been chosen Gaussian. If the simplest Gaussian is chosen, as in 
[4], it is of form 

P(X) = P(\X\) = 2n X p( X ) = %e-* 2 '* 2 (El) 

which is normalized and has the correct second moment (x 2 ) = a? The ambiguity is more 
obvious for the K model. Requiring that p(x) should be a monotonic, bounded function, 
exhibit an asymptotic if-Bessel function behavior and, for convenience, be amenable to 
analytic analysis we can choose 

P(X) = A^P- (E2) 

where Z = ^(x 2 + P 2 ) 1 ^ 2 - For any v > 3/4, values of A and (5 can be determined analytically 
[20] such as to satisfy the normalization and second moment criteria. For v — 1: (3 = va 
and A = [iraKo^)}' 1 where v = 0.1657 is the solution to K (v) = 2vK 1 (v). For small 
values of a, the K model, as compared with the Gauss model, describes a surface with 
higher probability of large slope (x » ex) at the expense of areas with a very small slope. 

In the macroscopic-roughness model the probability of scattering (local reflection) from 
f2j to f2, i.e. from angles (#«,0) to (0,ip), at suitably oriented surface elements is : 

dP 

—dQ = p( X )d 2 X = p(x)dXxd Xy (E3) 
The mapping x — > ^ involves the local polar angle of incidence, 9', where 

d = cos & = (c\ cos 2 (I) + c 2 sin 2 (|) ) ^ (E4) 
with c± = cos[(# ± 00/2]. We obtain 

dQ = ^fd 2 X (E5) 
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Inserting (E5) and (El) into (E3) gives the scattering distribution. For the Gaussian model, 
t in n^ dP 1 c ' 2 

where the in-plane momentum transfer q is given in (11). 

Eq. (E6) is consistent with the limiting case analyzed in ref. [4] (Eq. (26) of [4]) and it 
closely resembles the micro-roughness result (14) for /m) (for r] — 0). It exhibits a similar 
width of scattering distribution around the regular reflection angle 0) but no exactly 
specular intensity. Local angle (#') dependent loss could be incorporated in the form Ho(9') 
as for a plane surface, but the result is not identical to the micro-roughness result in the 
limit a « 1. 

Expression (E6) and similar expressions for non-Gaussian models are symmetric in and 
Q, thus do not satisfy the detailed balance requirement discussed in Section 6, and indeed 
give rise to loss of isotropy in computer simulations. Moreover, (E6) is catastrophically 
inadequate for glancing-angle incidence (0j — > n/2) due to shadowing and multiple reflection 
within a rough surface as discussed in Section 6. Although this angular range is small, it 
makes a large contribution to the mean value for isotropic incident UCN flux and cannot be 
neglected. 



(2akc + C-Y 



(E6) 
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